package scu.maqiang.numeric;

import scu.maqiang.csparsej.tdouble.*;
import scu.maqiang.cxsparsej.*;

import java.util.Arrays;

/**
 * 稀疏矩阵直接求解器，包装调用了csparsej包中的直接求解器, 可实现实系数与复系数稀疏矩阵直接求解
 * @author 马强
 *
 */
public class DirectSCSolver extends Solver{
    DZcs_common.DZcs matrix = null;
    DZcs_common.DZcss S = null;
    DZcs_common.DZcsn N = null;
    int n = 0;

    public DirectSCSolver() {
    }

    /**
     * 构造函数, 指定求解器要操作的实稀疏矩阵
     * @param mat 实稀疏矩阵
     */
    public DirectSCSolver(SCMatrix mat) {
        DZcs_common.DZcs temp = mat.toDZcs_matrix();
        matrix = DZcs_compress.cs_compress(temp);
        n = mat.getM();
    }

    public void setMatrix(SCMatrix mat) {
        DZcs_common.DZcs temp = mat.toDZcs_matrix();
        matrix = DZcs_compress.cs_compress(temp);
        //DZcs_print.cs_print(matrix, false);
        //System.exit(1);
        n = mat.getM();
    }

    /**
     * 使用LU分解方法求解一般系数矩阵的线性方程组
     * @param RHS 线性方程组右端项
     * @param order 使用的排序方法 (0 to 3)
     * @param tol 部分选主元的容许误差
     * @return 方程组的解
     */
    public Complex[] LUSolve(Complex[] RHS, int order, double tol) {
        DZcs_common.DZcsa result = new DZcs_common.DZcsa(RHS.length);
        for(int i = 0; i < RHS.length; i++) {
            result.x[2 * i] = RHS[i].real;
            result.x[2 * i + 1] = RHS[i].image;
        }
        long begin = System.nanoTime();
        DZcs_lusol.cs_lusol(order, matrix, result, tol);
        long end = System.nanoTime();
        Complex[] sol = Complex.zeros(RHS.length);
        for(int i = 0; i < RHS.length; i++) {
            sol[i].real = result.x[2 * i];
            sol[i].image = result.x[2 * i + 1];
        }
        //System.out.printf("min: %20.10e\t max: %20.10e\n", MVO.min(sol), MVO.max(result));
        System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
        return sol;
    }

    /**
     * 对一般系数矩阵进行LU分解, 与Solve_LU一起使用
     * @param order 排序方法 (0 to 3)
     * @see #Solve_LU(double[])
     * @return 是否分解成功
     */
//    public boolean Factorize_LU(int order, double tol) {
//        int n;
//        boolean ok;
//        if (!DZcs_util.CS_CSC(matrix))
//            return (false); /* check inputs */
//        S = DZcs_sqr.cs_sqr(order, matrix, false); /* ordering and symbolic analysis */
//        N = DZcs_lu.cs_lu(matrix, S, tol); /* numeric LU factorization */
//        ok = (S != null && N != null);
//        //System.out.println(Arrays.toString(b));
//        return (ok);
//    }

//    /**
//     * 对已经LU分解好的系数矩阵进行求解, 与Factorize_LU一起使用
//     * @param b 线性方程组右端项
//     * @see #Factorize_LU(int, double)
//     * @return 方程组解
//     */
//    public double[] Solve_LU(double[] b, int verbosity) {
//        int n = matrix.n;
//        double[] x = new double[n]; /* get workspace */
//        double bb[] = b.clone();
//        boolean ok = (S != null && N != null && b != null);
//        long begin = System.nanoTime();
//        if (ok) {
//            DZcs_ipvec.cs_ipvec(N.pinv, bb, x, n); /* x = b(p) */
//            DZcs_lsolve.cs_lsolve(N.L, x); /* x = L\x */
//            DZcs_usolve.cs_usolve(N.U, x); /* x = U\x */
//            DZcs_ipvec.cs_ipvec(S.q, x, bb, n); /* b(q) = x */
//        }
//        long end = System.nanoTime();
//        if(verbosity >= 1) {
//            System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(bb), MVO.max(bb));
//            System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//        }
//        return bb;
//    }

//	/**
//	 * 当使用这个构造函数时，求解问题使用PPLUSolve
//	 * @param A
//	 */
//	public DirectSSolver(SRMatrix A) {
//		this.A = A;
//		int m = A.getM();
//		P = A.symamd();
//		L = new SRMatrix(m);
//		U = new SRMatrix(m);
//		//P.setDiag(1.0);
//		L.setDiag(1.0);
//		LU();
//	}
//
//	/**
//	 *  当使用这个构造函数时，求解问题使用PLUSolve
//	 * @param A
//	 */
//	public DirectSSolver(SRMatrix A, int nBand) {
//		this.A = A;
//		this.nBand = nBand;
//		int m = A.getM();
//		//P = new SRMatrix(m);
//		L = new SRMatrix(m);
//		U = new SRMatrix(m);
//		//P.setDiag(1.0);
//		L.setDiag(1.0);
//		P = new int[m];
//		Arrays.setAll(P, i -> i);
//		PLU();
//		//System.out.println(P);
//	}
//
//	/**
//	 * 对稀疏矩阵A进行LU分解
//	 */
//	public void LU() {
//		long begin = System.nanoTime();
//		int m = A.getM();
//		int fillinCount = 0;
//		for(int i = 0; i < m - 1; i++) {
//			double maxval = A.getElement(i, i);
//			for(int j = i + 1; j < m; j++) {
//				SREntry ji = A.getSREntry(j, i);
//				if (ji != null) {
//					double ratio = ji.val / maxval;
//					ji.val = ratio;
//					for(SREntry ele : A.mat[i]) {
//						if (ele.col > i) {
//							SREntry jCol = A.getSREntry(j, ele.col);
//							if(jCol != null) {
//								jCol.val -= ele.val * ratio;
//							} else {
//								A.setElementFast(j, ele.col, -ele.val * ratio);
//								System.out.println("Fill in in\t " + j + "\t, " + ele.col);
//								fillinCount++;
//							}
//						}
//					}
//				}
//			}
//		}
//		System.out.println("Factorizing finished....");
//		System.out.println("Fill in number: " + fillinCount);
//		long end = System.nanoTime();
//		System.out.println("Factorizing elapsed...\t" + (end - begin) / 1.0e9 + " s");
//		//System.out.println("P:" + Arrays.toString(P));
//		//System.out.println("P:" + MVO.toString(P.toDenseMatrix()));
////		System.out.println("A:" + MVO.toString(A.toDenseMatrix()));
//		for(int i = 0; i < m; i++) {
//			for(SREntry ele : A.mat[i]) {
//				if (ele.col >= i) {
//					U.mat[i].add(new SREntry(ele.col, ele.val));
//				} else {
//					L.mat[i].add(new SREntry(ele.col, ele.val));
//				}
//			}
//		}
//	}
//
//	public void PPLUSolve(double[] RHS, double[] x, int verbosity) {
//		long begin = System.nanoTime();
//		int m = L.getM();
//		double[] Pb = new double[m];
//		Arrays.setAll(Pb, i -> RHS[P[i]]);
//		double[] C = new double[m];
//		C[0] = Pb[0];
//		for(int i = 1; i < m; i++) {
//			double temp = 0.0;
//			for(SREntry ele : L.mat[i]) {
//				if (ele.col < i) {
//					temp += C[ele.col] * ele.val;
//				}
//			}
//			C[i] = Pb[i] - temp;
//		}
//
//		double[] xx = new double[m];
//		double vMM = U.getElement(m - 1, m - 1);
//		xx[m - 1] = C[m - 1] / vMM;
//		for(int i = m - 2; i >=0; i--) {
//			double temp = 0.0;
//			for(SREntry ele : U.mat[i]) {
//				if (ele.col > i) {
//					temp += xx[ele.col] * ele.val;
//				} else {
//					vMM = ele.val;
//				}
//			}
//			xx[i] = (C[i] - temp) / vMM;
//		}
//		for(int i = 0; i < m; i++) {
//			x[P[i]] = xx[i];
//		}
//
//		if (verbosity >= 1) {
//			System.out.printf("min: %20.10e\t max: %20.10e\n", MVO.min(x), MVO.max(x));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//	}
//
//
//
//	/**
//	 * 对稀疏矩阵A进行PLU分解
//	 */
//	public void PLU() {
//		long begin = System.nanoTime();
//		int m = A.getM();
//		int fillinCount = 0;
//		for(int i = 0; i < m - 1; i++) {
//			Pair<Integer, Double> maxIdxResult = A.maxAbsColValDown(i, nBand);
//			int idxJ = maxIdxResult.getFirst();
//			double maxval = maxIdxResult.getSecond();
//			if (Math.abs(maxval) < Constants.Er) {
//				throw new ArithmeticException("Zero pivot element!");
//			}
//			if (maxIdxResult.getFirst() != i) {
//				A.exchangeRow(i, idxJ);
//				//System.out.println("OK!");
//				//P.exchangeRow(i, idxJ);
//				int temp = P[i];
//				P[i] = P[idxJ];
//				P[idxJ] = temp;
//			}
//			int end = Math.min(i + nBand + 3, m);
//
//			for(int j = i + 1; j < end; j++) {
//				SREntry ji = A.getSREntry(j, i);
//				if (ji != null) {
//					double ratio = ji.val / maxval;
//					ji.val = ratio;
//					for(SREntry ele : A.mat[i]) {
//						if (ele.col > i) {
//							SREntry jCol = A.getSREntry(j, ele.col);
//							if(jCol != null) {
//								jCol.val -= ele.val * ratio;
//							} else {
//								A.setElementFast(j, ele.col, -ele.val * ratio);
//								fillinCount++;
//								System.out.println("Fill in in\t " + j + "\t, " + ele.col);
//							}
//						}
//					}
//				}
//			}
//		}
//		System.out.println("Factorizing finished....");
//		System.out.println("Fill in number: " + fillinCount);
//		long end = System.nanoTime();
//		System.out.println("Factorizing elapsed...\t" + (end - begin) / 1.0e9 + " s");
//		//System.out.println("P:" + Arrays.toString(P));
//		//System.out.println("P:" + MVO.toString(P.toDenseMatrix()));
////		System.out.println("A:" + MVO.toString(A.toDenseMatrix()));
//		for(int i = 0; i < m; i++) {
//			for(SREntry ele : A.mat[i]) {
//				if (ele.col >= i) {
//					U.mat[i].add(new SREntry(ele.col, ele.val));
//				} else {
//					L.mat[i].add(new SREntry(ele.col, ele.val));
//				}
//			}
//		}
//
////		System.out.println(Arrays.toString(P));
////		System.out.println("P:" + MVO.toString(P.toDenseMatrix()));
////		System.out.println("L:" + MVO.toString(L.toDenseMatrix()));
////		System.out.println("U:" + MVO.toString(U.toDenseMatrix()));
//	}
//
//	/**
//	 *
//	 * @param RHS
//	 * @param x
//	 * @param verbosity
//	 */
//	public void PLUSolve(double[] RHS, double[] x, int verbosity) {
//		long begin = System.nanoTime();
//		int m = L.getM();
//		double[] Pb = new double[m];
//		Arrays.setAll(Pb, i -> RHS[P[i]]);
//		double[] C = new double[m];
//		C[0] = Pb[0];
//		for(int i = 1; i < m; i++) {
//			double temp = 0.0;
//			for(SREntry ele : L.mat[i]) {
//				if (ele.col < i) {
//					temp += C[ele.col] * ele.val;
//				}
//			}
//			C[i] = Pb[i] - temp;
//		}
//
//		double vMM = U.getElement(m - 1, m - 1);
//		x[m - 1] = C[m - 1] / vMM;
//		for(int i = m - 2; i >=0; i--) {
//			double temp = 0.0;
//			for(SREntry ele : U.mat[i]) {
//				if (ele.col > i) {
//					temp += x[ele.col] * ele.val;
//				} else {
//					vMM = ele.val;
//				}
//			}
//			x[i] = (C[i] - temp) / vMM;
//		}
//
//		if (verbosity >= 1) {
//			System.out.printf("min: %20.10e\t max: %20.10e\n", MVO.min(x), MVO.max(x));
//			long end = System.nanoTime();
//			System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");
//			System.out.println();
//		}
//	}
//
//
//	SRMatrix A;
//	int[] P;
//	SRMatrix L;
//	SRMatrix U;
//	int nBand;
}
